import numpy as np
import matplotlib.pyplot as plt
import plumed
import os

plt.rcParams['figure.dpi'] = 200

home_dir=os.getcwd()
print(home_dir)

/Users/valsson/WorkTmp-4/masterclass-22-11/Solution

# kB*T at 300 K in kJ/mol
kBT=2.494353

Exercise 1 - Biasing with one CV

It is assumed that exercise has already been run.

Time evolution of the CVs

Distance CV

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],distance[::4],'-')
plt.xlabel("Time [ps]")
plt.ylabel("Na-Cl distance [nm]")
plt.ylim([0.2,0.7])
#plt.xlim([0,1000])

(0.2, 0.7)

As we can see, the biased distance CV is diffusive during the biased simulation

Solvation of the Na atom (coordination number)

# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],coord[::4],'-')
plt.xlabel("Time [ps]")
plt.ylabel("Solvation of Na")

Text(0, 0.5, 'Solvation of Na')

Coupling between distance and solvation CV

# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(distance[::4],coord[::4],'.')
plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Solvation of Na")

Text(0, 0.5, 'Solvation of Na')

Time evoluation of the coeffients

%%bash
# Extract the time series of the coefficents
cd Exercise-1/
for i in `seq 25`
do
  ./ExtractCoeff.sh ${i} > coeff.${i}.data
done

for i in range(1,26):
    coeff_data = np.loadtxt("./Exercise-1/coeff.{0}.data".format(i))
    plt.plot(coeff_data[:,1],label="Average Coeff")
    plt.plot(coeff_data[:,2],label="Instantaneous Coeff")
    plt.ylabel("Coefficent".format(i))
    plt.legend()
    plt.title("Coefficent {0}".format(i))
    plt.show()

Convergence of the FES

# Each iteration is 1 ps
for iter in range(500,4001,500):
    fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.title("Free Energy Surface")
    plt.legend()
    plt.ylim([0,22])

# Each iteration is 1 ps
for iter in range(3000,4001,100):
    fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

Convergence of free energy difference

last_iter=4000
fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
#plt.plot(fes_data[:,0],fes_data[:,1]/kBT)
plt.plot(fes_data[:,0],np.exp(-fes_data[:,1]/kBT))
plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Probability Distribution")
plt.title("Probability Distribution")
plt.ylim([0,1.1])
barrier_location=0.36
plt.axvspan(0.2,barrier_location,alpha=0.2,color='orange')
plt.axvspan(barrier_location,0.7,alpha=0.2,color='green')

<matplotlib.patches.Polygon at 0x13e601f60>

def calc_free_energy_difference(fes_file):
    barrier_location=0.36
    fes_data = np.loadtxt(fes_file)
    distance = fes_data[:,0]
    fes = fes_data[:,1]
    prob = np.exp(-fes/kBT)
    prob_A = 0.0
    prob_B = 0.0
    for i in range(prob.size):
        if(distance[i]<barrier_location): 
            prob_A += prob[i]
        if(distance[i]>barrier_location): 
            prob_B += prob[i]
    free_energy_difference = -np.log(prob_A/prob_B)
    return free_energy_difference

iteration_ex1=[]
free_energy_difference_ex1=[]

for iter in range(100,4001,100):
    iteration_ex1.append(iter)
    fes_filename="./Exercise-1/fes.ves.iter-{0}.data".format(iter)
    fe_diff = calc_free_energy_difference(fes_filename)
    free_energy_difference_ex1.append(fe_diff)

plt.plot(iteration_ex1,free_energy_difference_ex1)
plt.xlabel("Time [ps]")
plt.ylabel("Free Energy Difference [kBT]")
plt.title("Free Energy Difference")
plt.ylim([-2.0,2.0])

(-2.0, 2.0)

Convergence of the bias potential and target distribution

# Each iteration is 1 ps
for iter in range(1000,4001,1000):
    bias_data = np.loadtxt("./Exercise-1/bias.ves.iter-{0}.data".format(iter))
    plt.plot(bias_data[:,0],bias_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Bias [kBT]")
    plt.legend()
    plt.title("Bias Potential")
    # plt.ylim([0,22])

# Each iteration is 1 ps
for iter in range(1000,4001,1000):
    bias_data = np.loadtxt("./Exercise-1/targetdist.ves.iter-{0}.data".format(iter))
    plt.plot(bias_data[:,0],bias_data[:,1],label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Target Distribution")
    plt.legend()
    plt.title("Target Distribution")

# Each iteration is 1 ps
for iter in range(1000,4001,1000):
    bias_data = np.loadtxt("./Exercise-1/targetdist.ves.log.iter-{0}.data".format(iter))
    plt.plot(bias_data[:,0],bias_data[:,1],label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Negative Logarithm of Target Distribution [kBT]")
    plt.legend()
    plt.title("Negative Logarithm of Target Distribution")
    plt.ylim([0,2.5])

The negative logarithm of the target distribution is the effective that we are sampling on. As we are using a well-tempered target distribution with a bias factor of 10, this should be an effective FES where the barriers have been reduced by a factor of 10. As we can see in the FES figure above, the barrier is around 5 kBT, while in this figure the effective barriers are around 0.5 kBT so this matches.

Exercise 2 - Reweighthing

Behaviour of bias potential and weights

First we should look at the time series of the bias potential and corresponding weights to check how much we should trim from the beginning of the colvar file

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
bias = colvar_data[:,3]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],bias[::4],'-')
plt.xlabel("Time [ps]")
plt.ylabel("Bias [kJ/mol]")

Text(0, 0.5, 'Bias [kJ/mol]')

weight = np.exp(bias-np.max(bias))
plt.plot(time[::1],weight[::1],'-')
plt.xlabel("Time [ps]")
plt.ylabel("Weight")

Text(0, 0.5, 'Weight')

plt.plot(distance[::1],weight[::1],'.')
plt.xlabel("Na-Cl distance")
plt.ylabel("Weight")

Text(0, 0.5, 'Weight')

Based on this, we should trim the first 500-1000 ps from the colvar file

Reweighthing on Na-Cl distance

%%bash
cat ./Exercise-2/plumed_reweight.dat

# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist
  GRID_MIN=0.2
  GRID_MAX=0.9
  GRID_BIN=60
  KERNEL=DISCRETE
  LOGWEIGHTS=weights
  LABEL=hg_dist
... HISTOGRAM

fes_dist: CONVERT_TO_FES GRID=hg_dist TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist FILE=fes-reweight.dist.data FMT=%24.16e

os.chdir(home_dir)
os.chdir("./Exercise-2")
for i in range(0,1001,100):
    os.system("./trim-colvar-file.py --colvar-file ../Exercise-1/colvar.data --output-file colvar_reweight.data --time-min {0} &> /dev/null ".format(i))
    os.system("plumed driver --plumed plumed_reweight.dat --noatoms &> /dev/null ")
    fes_data=np.loadtxt("fes-reweight.dist.data")
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Trim {} ps".format(i))
    plt.xlabel("Na-Cl distance [nm]")  
    plt.ylabel("Free Energy [kBT]")
    plt.title("Reweighted FES")
    plt.legend()
    os.system("rm -f fes-reweight.dist.data")
    os.system("rm -f colvar_reweight.data")
    plt.ylim([0,18])
os.chdir(home_dir)

os.chdir(home_dir)
os.chdir("./Exercise-2")
os.system("rm -f bck*.fes-reweight.dist.data")
os.system("rm -f fes-reweight.dist.data")
trim_time="500"
os.system("./trim-colvar-file.py --colvar-file ../Exercise-1/colvar.data --output-file colvar_reweight.data --time-min {0} &> /dev/null ".format(trim_time))
os.system("plumed driver --plumed plumed_reweight.dat --noatoms &> /dev/null ")
os.chdir(home_dir)

Comparing FES directly from bias potential and from reweighting

last_iter=4000
fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="FES from bias")

fes_rew_data = np.loadtxt("./Exercise-2/fes-reweight.dist.data")
plt.plot(fes_rew_data[:,0],fes_rew_data[:,1]/kBT,label="FES from reweighting ")

plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from bias and reweighting")
plt.ylim([0,15])

(0.0, 15.0)

Smoothing using kernel density estimation

%%bash
cd Exercise-2/
cat << EOF > plumed_reweight.kde.dat 
# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist
  GRID_MIN=0.2
  GRID_MAX=0.9
  GRID_BIN=200
  BANDWIDTH=0.004
  LOGWEIGHTS=weights
  LABEL=hg_dist
... HISTOGRAM

fes_dist: CONVERT_TO_FES GRID=hg_dist TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist FILE=fes-reweight.kde.dist.data FMT=%24.16e
EOF

plumed driver --plumed plumed_reweight.kde.dat --noatoms &> /dev/null
rm -f bck.*.fes-reweight.kde.dist.data
ls 

cd ..

colvar_reweight.data
fes-reweight.dist-coord.data
fes-reweight.dist.data
fes-reweight.kde.dist.data
plumed_reweight.2D.dat
plumed_reweight.dat
plumed_reweight.kde.dat
trim-colvar-file.py

last_iter=4000
fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="FES from bias")

fes_rew_data = np.loadtxt("./Exercise-2/fes-reweight.dist.data")
plt.plot(fes_rew_data[:,0],fes_rew_data[:,1]/kBT,label="FES from reweighting ")

fes_rew_kde_data = np.loadtxt("./Exercise-2/fes-reweight.kde.dist.data")
plt.plot(fes_rew_kde_data[:,0],fes_rew_kde_data[:,1]/kBT,label="FES from reweighting - KDE")

plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from bias and reweighting")
plt.ylim([0,15])

(0.0, 15.0)

Reweighthing on Na-Cl distance and solvation CV

%%bash
cd Exercise-2/
cat << EOF > plumed_reweight.2D.dat 
# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist,coord
  GRID_MIN=0.2,2.5
  GRID_MAX=0.9,7.5
  GRID_BIN=200,200
  BANDWIDTH=0.004,0.04
  LOGWEIGHTS=weights
  LABEL=hg_dist_coord
... HISTOGRAM

fes_dist_coord: CONVERT_TO_FES GRID=hg_dist_coord TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist_coord FILE=fes-reweight.dist-coord.data FMT=%24.16e
EOF

plumed driver --plumed plumed_reweight.2D.dat --noatoms &> /dev/null
rm -f bck.*.fes-reweight.dist-coord.data
ls 

cd ..

colvar_reweight.data
fes-reweight.dist-coord.data
fes-reweight.dist.data
fes-reweight.kde.dist.data
plumed_reweight.2D.dat
plumed_reweight.dat
plumed_reweight.kde.dat
trim-colvar-file.py

fes_data = np.loadtxt("./Exercise-2/fes-reweight.dist-coord.data")
distance = fes_data[:,0].reshape(201,201)
coord =    fes_data[:,1].reshape(201,201)
fes_2d =   fes_data[:,2].reshape(201,201)
from matplotlib import cm
plt.contourf(distance,coord,fes_2d/kBT, levels=np.linspace(0,10,101), cmap=cm.jet)
plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Solvation of Na")
plt.title("FES - Reweighting")
plt.colorbar(label="Free Energy [kBT]", ticks=range(0,11))

<matplotlib.colorbar.Colorbar at 0x13eb4ffa0>

Exercise 3 - Another independent simulation (biasing the distance)

colvar_data = np.loadtxt("./Exercise-3/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],distance[::4],'-')
plt.xlabel("Time [ps]")
plt.ylabel("Na-Cl distance [nm]")
plt.ylim([0.2,0.7])

(0.2, 0.7)

colvar_data = np.loadtxt("./Exercise-3/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],distance[::4],'-',label="Run 1")
plt.xlabel("Time [ps]")
plt.ylabel("Na-Cl distance [nm]")
plt.ylim([0.2,0.7])

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(time[::4],distance[::4],'--',label="Run 2")
plt.legend()

<matplotlib.legend.Legend at 0x13ec07190>

# Each iteration is 1 ps
for iter in range(500,4001,500):
    fes_data = np.loadtxt("./Exercise-3/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

# Each iteration is 1 ps
for iter in range(3000,4001,100):
    fes_data = np.loadtxt("./Exercise-3/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

last_iter=4000

fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1")

fes_data = np.loadtxt("./Exercise-3/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 2")


plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from two different independent runs")
plt.ylim([0,15])

(0.0, 15.0)

%%bash
cd Exercise-3/
cat << EOF > plumed_reweight.kde.dat 
# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist
  GRID_MIN=0.2
  GRID_MAX=0.9
  GRID_BIN=200
  BANDWIDTH=0.004
  LOGWEIGHTS=weights
  LABEL=hg_dist
... HISTOGRAM

fes_dist: CONVERT_TO_FES GRID=hg_dist TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist FILE=fes-reweight.kde.dist.data FMT=%24.16e
EOF

trim_time=500 
../Exercise-2/trim-colvar-file.py --colvar-file ./colvar.data --output-file colvar_reweight.data --time-min ${trim_time} &> /dev/null

plumed driver --plumed plumed_reweight.kde.dat --noatoms &> /dev/null
rm -f bck.*.fes-reweight.kde.dist.data


cd ..

last_iter=4000


fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Bias")

fes_data = np.loadtxt("./Exercise-2/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Reweight",linestyle='--')

fes_data = np.loadtxt("./Exercise-3/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 2 - Bias")

fes_data = np.loadtxt("./Exercise-3/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 2 - Reweight",linestyle='--')

plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from two different independent runs")
plt.ylim([0,15])

(0.0, 15.0)

iteration_ex3=[]
free_energy_difference_ex3=[]

for iter in range(100,4001,100):
    iteration_ex3.append(iter)
    fes_filename="./Exercise-3/fes.ves.iter-{0}.data".format(iter)
    fe_diff = calc_free_energy_difference(fes_filename)
    free_energy_difference_ex3.append(fe_diff)

plt.plot(iteration_ex1,free_energy_difference_ex1,label="Run 1 (Ex 1)")
plt.plot(iteration_ex3,free_energy_difference_ex3,label="Run 2 (Ex 3)")
plt.xlabel("Time [ps]")
plt.ylabel("Free Energy Difference [kBT]")
plt.title("Free Energy Difference")
plt.ylim([-2.0,2.0])
plt.legend()

<matplotlib.legend.Legend at 0x13e7e0940>

Exercise 4 - Biasing with two CVs

colvar_data = np.loadtxt("./Exercise-4/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],distance[::4],'-')
plt.xlabel("Time [ps]")
plt.ylabel("Na-Cl distance [nm]")
plt.ylim([0.2,0.7])

(0.2, 0.7)

plt.plot(time[::4],coord[::4],'-')
plt.xlabel("Time [ps]")
plt.ylabel("solvation of Na")

Text(0, 0.5, 'solvation of Na')

colvar_data = np.loadtxt("./Exercise-4/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(time[::4],coord[::4],'-',label="Biasing distance/solvation")

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(time[::4],coord[::4],'-',label="Biasing distance")

plt.xlabel("Time [ps]")
plt.ylabel("solvation of Na")
plt.legend()

<matplotlib.legend.Legend at 0x13ec6f340>

colvar_data = np.loadtxt("./Exercise-4/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(distance[::4],coord[::4],'.',label="Biasing distance/solvation")

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(distance[::4],coord[::4],'.',label="Biasing distance")

plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("solvation of Na")
plt.legend()

<matplotlib.legend.Legend at 0x13ed93c40>

# Each iteration is 1 ps
for iter in range(500,4001,500):
    fes_data = np.loadtxt("./Exercise-4/fes.ves.proj-1.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

# Each iteration is 1 ps
for iter in range(3000,4001,100):
    fes_data = np.loadtxt("./Exercise-4/fes.ves.proj-1.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

# Each iteration is 1 ps
for iter in range(3000,4001,100):
    fes_data = np.loadtxt("./Exercise-4/fes.ves.proj-2.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Solvation of Na")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,6])

last_iter=4000

fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Biasing distance")

fes_data = np.loadtxt("./Exercise-3/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 2 - Biasing distance")

fes_data = np.loadtxt("./Exercise-4/fes.ves.proj-1.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 3 - Biasing distance/solvation")


plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from different runs")
plt.ylim([0,10])

(0.0, 10.0)

%%bash
cd Exercise-4/
cat << EOF > plumed_reweight.kde.dat 
# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist
  GRID_MIN=0.2
  GRID_MAX=0.9
  GRID_BIN=200
  BANDWIDTH=0.004
  LOGWEIGHTS=weights
  LABEL=hg_dist
... HISTOGRAM

fes_dist: CONVERT_TO_FES GRID=hg_dist TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist FILE=fes-reweight.kde.dist.data FMT=%24.16e
EOF

trim_time=500 
../Exercise-2/trim-colvar-file.py --colvar-file ./colvar.data --output-file colvar_reweight.data --time-min ${trim_time} &> /dev/null

plumed driver --plumed plumed_reweight.kde.dat --noatoms &> /dev/null
rm -f bck.*.fes-reweight.kde.dist.data


cd ..

last_iter=4000

fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Biasing distance")

fes_data = np.loadtxt("./Exercise-2/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Reweight",linestyle='--')

fes_data = np.loadtxt("./Exercise-3/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 2 - Biasing distance")

fes_data = np.loadtxt("./Exercise-3/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 2 - Reweight",linestyle='--')

fes_data = np.loadtxt("./Exercise-4/fes.ves.proj-1.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 3 - Biasing distance/solvation")

fes_data = np.loadtxt("./Exercise-4/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 3 - Reweight",linestyle='--')

plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from different runs")
plt.ylim([0,10])

(0.0, 10.0)

iteration_ex4=[]
free_energy_difference_ex4=[]

for iter in range(100,4001,100):
    iteration_ex4.append(iter)
    fes_filename="./Exercise-4/fes.ves.proj-1.iter-{0}.data".format(iter)
    fe_diff = calc_free_energy_difference(fes_filename)
    free_energy_difference_ex4.append(fe_diff)

plt.plot(iteration_ex1,free_energy_difference_ex1,label="Run 1 (Ex 1)")
plt.plot(iteration_ex3,free_energy_difference_ex3,label="Run 2 (Ex 3)")
plt.plot(iteration_ex4,free_energy_difference_ex4,label="Run 3 (Ex 4)")
plt.xlabel("Time [ps]")
plt.ylabel("Free Energy Difference [kBT]")
plt.title("Free Energy Difference")
plt.ylim([-2.0,2.0])
plt.legend()

<matplotlib.legend.Legend at 0x13e8d3790>

last_iter=4000
fes_data = np.loadtxt("./Exercise-4/fes.ves.iter-{0}.data".format(last_iter))
distance = fes_data[:,0].reshape(301,301)
coord =    fes_data[:,1].reshape(301,301)
fes_2d =   fes_data[:,2].reshape(301,301)
from matplotlib import cm
plt.contourf(distance,coord,fes_2d/kBT, levels=np.linspace(0,7,71), cmap=cm.jet)
plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Solvation of Na")
plt.title("FES - Biasing distance and solvation CVs")
plt.colorbar(label="Free Energy [kBT]", ticks=range(0,8))

<matplotlib.colorbar.Colorbar at 0x13ef59570>

for iter in range(1000,4001,1000):
    fes_data = np.loadtxt("./Exercise-4/fes.ves.iter-{0}.data".format(iter))
    distance = fes_data[:,0].reshape(301,301)
    coord =    fes_data[:,1].reshape(301,301)
    fes_2d =   fes_data[:,2].reshape(301,301)
    from matplotlib import cm
    plt.contourf(distance,coord,fes_2d/kBT, levels=np.linspace(0,7,71), cmap=cm.jet)
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Solvation of Na")
    plt.title("FES - Iteration {0}".format(iter))
    plt.colorbar(label="Free Energy [kBT]", ticks=range(0,8))
    plt.show()

Exercise 6 - Uniform target distribution

colvar_data = np.loadtxt("./Exercise-6/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],distance[::4],'-',label="Uniform target distribution")
plt.xlabel("Time [ps]")
plt.ylabel("Na-Cl distance [nm]")
plt.ylim([0.2,0.7])

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(time[::4],distance[::4],'--',label="Well-tempered target distribution")
plt.legend()

<matplotlib.legend.Legend at 0x13e36c970>

# Each iteration is 1 ps
for iter in range(500,4001,500):
    fes_data = np.loadtxt("./Exercise-6/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

# Each iteration is 1 ps
for iter in range(3000,4001,100):
    fes_data = np.loadtxt("./Exercise-6/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

%%bash
cd Exercise-6/
cat << EOF > plumed_reweight.kde.dat 
# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist
  GRID_MIN=0.2
  GRID_MAX=0.9
  GRID_BIN=200
  BANDWIDTH=0.004
  LOGWEIGHTS=weights
  LABEL=hg_dist
... HISTOGRAM

fes_dist: CONVERT_TO_FES GRID=hg_dist TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist FILE=fes-reweight.kde.dist.data FMT=%24.16e
EOF

trim_time=500 
../Exercise-2/trim-colvar-file.py --colvar-file ./colvar.data --output-file colvar_reweight.data --time-min ${trim_time} &> /dev/null

plumed driver --plumed plumed_reweight.kde.dat --noatoms &> /dev/null
rm -f bck.*.fes-reweight.kde.dist.data


cd ..

last_iter=4000

fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Well-Tempered")

fes_data = np.loadtxt("./Exercise-2/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Reweight",linestyle='--')

fes_data = np.loadtxt("./Exercise-6/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 4 - Uniform")

fes_data = np.loadtxt("./Exercise-6/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 4 - Reweight",linestyle='--')



plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from different runs")
plt.ylim([0,10])

(0.0, 10.0)

iteration_ex6=[]
free_energy_difference_ex6=[]

for iter in range(100,4001,100):
    iteration_ex6.append(iter)
    fes_filename="./Exercise-6/fes.ves.iter-{0}.data".format(iter)
    fe_diff = calc_free_energy_difference(fes_filename)
    free_energy_difference_ex6.append(fe_diff)

plt.plot(iteration_ex1,free_energy_difference_ex1,label="Run 1 (Ex 1)")
plt.plot(iteration_ex3,free_energy_difference_ex3,label="Run 2 (Ex 3)")
plt.plot(iteration_ex4,free_energy_difference_ex4,label="Run 3 (Ex 4)")
plt.plot(iteration_ex6,free_energy_difference_ex6,label="Run 4 (Ex 6)")
plt.xlabel("Time [ps]")
plt.ylabel("Free Energy Difference [kBT]")
plt.title("Free Energy Difference")
plt.ylim([-2.0,2.0])
plt.legend()

<matplotlib.legend.Legend at 0x13f366ec0>

Exercise 7 - Legendre Polynomials

colvar_data = np.loadtxt("./Exercise-7/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
# the colvar file is outputted every 0.25 ps, only plot every 1.0 ps
plt.plot(time[::4],distance[::4],'-',label="Legendre Basis Functions")
plt.xlabel("Time [ps]")
plt.ylabel("Na-Cl distance [nm]")
plt.ylim([0.2,0.7])

colvar_data = np.loadtxt("./Exercise-1/colvar.data")
time = colvar_data[:,0]
distance = colvar_data[:,1]
coord = colvar_data[:,2]
plt.plot(time[::4],distance[::4],'--',label="Wavelet Basis Functions")
plt.legend()

<matplotlib.legend.Legend at 0x13f3100d0>

# Each iteration is 1 ps
for iter in range(500,4001,500):
    fes_data = np.loadtxt("./Exercise-7/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

# Each iteration is 1 ps
for iter in range(3000,4001,100):
    fes_data = np.loadtxt("./Exercise-7/fes.ves.iter-{0}.data".format(iter))
    plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="{0} ps".format(iter))
    plt.xlabel("Na-Cl distance [nm]")
    plt.ylabel("Free Energy [kBT]")
    plt.legend()
    plt.title("Free Energy Surface")
    plt.ylim([0,10])

%%bash
cd Exercise-7/
cat << EOF > plumed_reweight.kde.dat 
# vim:ft=plumed

dist:   READ FILE=colvar_reweight.data IGNORE_TIME VALUES=dist
coord:  READ FILE=colvar_reweight.data IGNORE_TIME VALUES=coord
ves:    READ FILE=colvar_reweight.data IGNORE_TIME VALUES=ves.bias

weights: REWEIGHT_BIAS TEMP=300 ARG=ves.bias

HISTOGRAM ...
  ARG=dist
  GRID_MIN=0.2
  GRID_MAX=0.9
  GRID_BIN=200
  BANDWIDTH=0.004
  LOGWEIGHTS=weights
  LABEL=hg_dist
... HISTOGRAM

fes_dist: CONVERT_TO_FES GRID=hg_dist TEMP=300 MINTOZERO
DUMPGRID GRID=fes_dist FILE=fes-reweight.kde.dist.data FMT=%24.16e
EOF

trim_time=500 
../Exercise-2/trim-colvar-file.py --colvar-file ./colvar.data --output-file colvar_reweight.data --time-min ${trim_time} &> /dev/null

plumed driver --plumed plumed_reweight.kde.dat --noatoms &> /dev/null
rm -f bck.*.fes-reweight.kde.dist.data


cd ..

last_iter=4000

fes_data = np.loadtxt("./Exercise-1/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Wavelet BF")

fes_data = np.loadtxt("./Exercise-2/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 1 - Reweight",linestyle='--')

fes_data = np.loadtxt("./Exercise-7/fes.ves.iter-{0}.data".format(last_iter))
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 5 - Legendre BF")

fes_data = np.loadtxt("./Exercise-7/fes-reweight.kde.dist.data")
plt.plot(fes_data[:,0],fes_data[:,1]/kBT,label="Run 5 - Reweight",linestyle='--')



plt.xlabel("Na-Cl distance [nm]")
plt.ylabel("Free Energy [kBT]")
plt.legend()
plt.title("FES from different runs")
plt.ylim([0,10])

(0.0, 10.0)

iteration_ex7=[]
free_energy_difference_ex7=[]

for iter in range(100,4001,100):
    iteration_ex7.append(iter)
    fes_filename="./Exercise-7/fes.ves.iter-{0}.data".format(iter)
    fe_diff = calc_free_energy_difference(fes_filename)
    free_energy_difference_ex7.append(fe_diff)

plt.plot(iteration_ex1,free_energy_difference_ex1,label="Run 1 (Ex 1)")
plt.plot(iteration_ex3,free_energy_difference_ex3,label="Run 2 (Ex 3)")
plt.plot(iteration_ex4,free_energy_difference_ex4,label="Run 3 (Ex 4)")
plt.plot(iteration_ex6,free_energy_difference_ex6,label="Run 4 (Ex 6)")
plt.plot(iteration_ex7,free_energy_difference_ex7,label="Run 5 (Ex 7)")
plt.xlabel("Time [ps]")
plt.ylabel("Free Energy Difference [kBT]")
plt.title("Free Energy Difference")
plt.ylim([-2.0,2.0])
plt.legend()

<matplotlib.legend.Legend at 0x13ee67700>